#pragma once
#include <vector>
#include <lapacke.h>
#include <iostream>
using namespace std;
class MatrixV
{
private:
    /// 几个参数的关系为 h = 1/(m+1), r = (k*nu)/(h*h)
    double r,h,k,nu;
    int m;
public:
    MatrixV(double _r, double _nu, double _h) : r(_r), nu(_nu), h(_h)
	{
	    m = (int)1/h - 1;
 	    k = h*h*r/nu;
	};
    MatrixV(double _nu, double _h) : nu(_nu), h(_h)
	{
	    m = (int)1/h - 1;
 	    r = 1;
	    k = 1;
	};
    vector<vector<double>> A()
	{
	    int i;
	    vector<vector<double>> A(m);
	    A[0].resize(m);
	    A[0][0] = -2*r/k;
	    for (i = 1; i < m; i++)
	    {
		A[i].resize(m);
		A[i][i] = A[0][0];
		A[i][i-1] = -A[0][0]/2;
		A[i-1][i] = -A[0][0]/2;
	    }
	    return A;
	};
    vector<vector<double>> A0()
	{
	    int i;
	    vector<vector<double>> A(9*(m+1)+m);
	    A[0].resize(9*(m+1)+m);
	    for (i = 1; i < 9*(m+1)+m; i++)
	    {
		A[i].resize(9*(m+1)+m);
		A[i][i-1] = -1;
		A[i-1][i] = 1;
	    }
	    return A;
	};
    vector<vector<double>> B0()
	{
	    int i;
	    vector<vector<double>> A(9*(m+1)+m);
	    A[0].resize(9*(m+1)+m);
	    for (i = 1; i < 9*(m+1)+m; i++)
	    {
		A[i].resize(9*(m+1)+m);
		A[i][i-1] = 1;
		A[i-1][i] = 1;
	    }
	    return A;
	};
    vector<vector<double>> C0()
	{
	    int i;
	    vector<vector<double>> A(9*(m+1)+m);
	    A[0].resize(9*(m+1)+m);
	    A[0][0] = -2;
	    for (i = 1; i < 9*(m+1)+m; i++)
	    {
		A[i].resize(9*(m+1)+m);
		A[i][i-1] = 1;
		A[i][i] = -2;
		A[i-1][i] = 1;
	    }
	    return A;
	};
    vector<vector<double>> E0()
	{
	    int i;
	    vector<vector<double>> A(9*(m+1)+m);
	    A[0].resize(9*(m+1)+m);
	    A[1].resize(9*(m+1)+m);
	    A[0][0] = 3;
	    A[1][0] = -4;
	    A[1][1] = 3;
	    for (i = 2; i < 9*(m+1)+m; i++)
	    {
		A[i].resize(9*(m+1)+m);
		A[i][i-1] = -4;
		A[i][i] = 3;
		A[i][i-2] = 1;
	    }
	    return A;
	};
    vector<vector<double>> F0()
	{
	    int i;
	    vector<vector<double>> A(9*(m+1)+m);
	    A[0].resize(9*(m+1)+m);
	    A[1].resize(9*(m+1)+m);
	    A[0][0] = 1;
	    A[1][0] = -2;
	    A[1][1] = 1;
	    for (i = 2; i < 9*(m+1)+m; i++)
	    {
		A[i].resize(9*(m+1)+m);
		A[i][i-1] = -2;
		A[i][i] = 1;
		A[i][i-2] = 1;
	    }
	    return A;
	};
    // vector<vector<double>> A0()
    // 	{
    // 	    int i;
    // 	    vector<vector<double>> A(m);
    // 	    A[0].resize(m);
    // 	    for (i = 1; i < m; i++)
    // 	    {
    // 		A[i].resize(m);
    // 		A[i][i-1] = -1;
    // 		A[i-1][i] = 1;
    // 	    }
    // 	    return A;
    // 	};
    vector<vector<double>> I()
	{
	    vector<vector<double>> I(m);
	    for (int i = 0; i < m; i++)
	    {
		I[i].resize(m);
		I[i][i] = 1;
	    }
	    return I;
	};
    vector<vector<double>> I0()
	{
	    vector<vector<double>> I(9*(m+1)+m);
	    for (int i = 0; i < 9*(m+1)+m; i++)
	    {
		I[i].resize(9*(m+1)+m);
		I[i][i] = 1;
	    }
	    return I;
	};
};
/// 向量加法
template <typename T>
vector<T> operator+(vector<T> V1, vector<T> V2)
{
    int n1 = V1.size(), n2 = V2.size();
    if (n1 != n2)
    {
 	cout << "向量加法的元素个数不相等" << endl;
 	exit(0);
    }
    for (int i = 0; i < n1; i++)
    {
	V1[i] = V1[i] + V2[i];
    }
    return V1;
};
/// 向量减法
template <typename T>
vector<T> operator-(vector<T> V1, vector<T> V2)
{
    int n1 = V1.size(), n2 = V2.size();
    if (n1 != n2)
    {
	cout << "向量减法的元素个数不相等" << endl;
	exit(0);
    }
    for (int i = 0; i < n1; i++)
    {
 	V1[i] = V1[i] - V2[i];
    }
    return V1;
};
/// 矩阵加法
template <typename T>
vector<vector<T>> operator+(vector<vector<T>> M1, vector<vector<T>> M2)
{
    int r1 = M1.size(), r2 = M2.size(), c1 = M1[0].size(), c2 = M2[0].size();
    if (r1 != r2 || c1 != c2)
    {
 	cout << "矩阵加法的行列不相等" << endl;
 	exit(0);
    }
    for (int i = 0; i < r1; i++)
    {
	M1[i] = M1[i] + M2[i];
    }
    return M1;
};
/// 矩阵减法
template <typename T>
vector<vector<T>> operator-(vector<vector<T>> M1, vector<vector<T>> M2)
{
    int r1 = M1.size(), r2 = M2.size(), c1 = M1[0].size(), c2 = M2[0].size();
    if (r1 != r2 || c1 != c2)
    {
 	cout << "矩阵减法的行列不相等" << endl;
 	exit(0);
    }
    for (int i = 0; i < r1; i++)
    {
	M1[i] = M1[i] - M2[i];
    }
    return M1;
};
/// 常数乘向量
template <typename T>
vector<T> operator*(double C, vector<T> V)
{
    int n = V.size();
    for (int i = 0; i < n; i++)
    {
 	V[i] = C*V[i];
    }
    return V;
};
/// 常数乘矩阵
template <typename T>
vector<vector<T>> operator*(double C, vector<vector<T>> M)
{
    int r = M.size();
    for (int i = 0; i < r; i++)
    {
 	M[i] = C*M[i];
    }
    return M;
};
/// 向量输出
template <typename T>
void print0(vector<T> V)
{
    int n = V.size();
    for (int i = 0; i < n; i++)
    {
	cout << V[i] << " ";
    }
    cout << ";" << endl;
};
/// 向量输出
template <typename T>
void print(vector<T> V)
{
    int n = V.size();
    cout << "[0, ";
    for (int i = 0; i < n; i++)
    {
	cout << V[i] << ", ";
    }
    cout << "0];" << endl;
};
/// 矩阵输出
template <typename T>
void print0(vector<vector<T>> M)
{
    int r = M.size(), c = M[0].size();
    for (int i = 0; i < r; i++)
    {
	print0(M[i]);
    }
};
/// 矩阵乘以向量
template <typename T>
vector<T> operator*(vector<vector<T>> M, vector<T> V)
{
    int r = M.size(), c = M[0].size(), n = V.size();
    vector<T> RV(r);
    if (c != n)
    {
	cout << "矩阵列数与向量长度不一致" << endl;
	exit(0);
    }
    for (int i = 0; i < r; i++)
    {
	for (int j = 0; j < n; j++)
	{
	    RV[i] = RV[i] + M[i][j]*V[j];
	}
    }
    return RV;
};
/// AU = F
vector<double> solveU(vector<vector<double>> A0, vector<double> F0)
{
    int m = A0.size(), n = A0[0].size(), k = F0.size();
    if (m != k)
    {
	cout << "length is not catch." << endl;
	exit(0);
    }
    vector<double> U(m);
    double *A = new double[m*m];
    double *F = new double[m];
    for (int j = 0; j < m; j++)
    {
	F[j] = F0[j];
	for (int l = 0; l < m; l++)
	{
	    A[j+l*m] = A0[j][l];
	}
    }
    lapack_int info,m0,n0,lda,ldb,nrhs;
    m0 = m;
    n0 = m;
    nrhs = 1;
    lda = m;
    ldb = m;
    info = LAPACKE_dgels(LAPACK_COL_MAJOR,'N',m0,n0,nrhs,A,lda,F,ldb);
    for(int i = 0; i < m; i++)
    {
	U[i] = F[i];
    }
    return U;
};

